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Abstract 

Polynomial systems occur in many fields of science and engineering. Polynomial homotopy 
continuation methods apply symbolic-numeric algorithms to solve polynomial systems. We 
describe the design and implementation of our web interface and reflect on the application of 
polynomial homotopy continuation methods to solve polynomial systems in the cloud. Via 
the graph isomorphism problem we organize and classify the polynomial systems we solved. 
The classification with the canonical form of a graph identifies newly submitted systems with 
systems that have already been solved. 
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1 Introduction 

The widespread availability and use of high speed internet connections combined with relatively 
inexpensive hardware enabled cloud computing. In cloud computing, users of software no longer 
download and install software, but connect via a browser to a web site, and interact with the 
software through a web interface. Computations happen at some remote server and the data 
(input as well as output) are stored and maintained remotely. Quoting [27], “Large, virtualized 
pools of computational resources raise the possibility of a new, advantageous computing paradigm 
for scientific research.” 

This model of computing offers several advantages to the user; we briefly mention three 
benefits. First, installing software can be complicated and a waste of time, especially if one wants 
to perform only one single experiment to check whether the software will do what is desired 
- the first author of [28j has an account on our web server. In cloud computing, the software 
installation is replaced with a simple sign up, as common as logging into a web store interface. 
One should not have to worry about upgrading installed software to newer versions. The second 
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advantage is that for computationally intensive tasks, the web server can be aided by a farm of 
compute servers. Thirdly, the input and output files are managed at the server. The user should 
not worry about storage, as the web server could be aided by file servers. A good web interface 
helps to formulate the input problems and manage the computed results. 

In this paper, we describe a prototype of a first web interface to the blackbox solver of PHC- 
pack m- This solver applies homotopy continuation methods to polynomial systems. Its black¬ 
box solver, available as phc -b at the command line, seems the most widely used and popular fea¬ 
ture of the software. Our web interface is currently running at https://kepler.math.uic.edu. 

In addition to the technical aspects of designing a web interface, we investigate what it 
means to run a blackbox solver in the cloud. Because the computations happen remotely, the 
blackbox solver not only hides the complexity of the algorithms, but also the actual cost of the 
computations. In principle, the solver could use just one single core of a traditional processor, 
or a distributed cluster of computers, accelerated with graphics processing units [31], [35]. One 
should consider classifying the hardness of an input problem and allocating the proper resources 
to solve the problem. This classification problem could be aided by mining a database of solved 
problems. 

To solve the classification problem we show that the problem of deciding whether two sets of 
support sets are isomorphic can be reduced to the graph isomorphism problem [21], [22]. Our 
classification problem is related to the isomorphism of polynomials problem [24] in multivariate 
cryptology 0, CD]. Support sets of polynomials span Newton polytopes and a related problems 
is the polytope isomorphism problem, see |18j for its complexity, which is as hard as the graph 
isomorphism problem. 


2 Related Work and Alternative Approaches 

The Sage notebook interface and the newer SageMathCloud can be alternative solutions to setting 
up a standalone cloud service. PHCpack is an optional package in Sage [30], available through the 
phc.py interface, developed by Marshall Hampton and Alex Jokela, based on the earlier efforts 
of Kathy Piret [26] and William Stein. We plan to upgrade the existing phc. py in Sage with 
phcpy [32]. 

The computational algebraic geometry software Macaulay2 [13] distributes PHCpack.m2 [14] . 
which is a package that interfaces to PHCpack. Macaulay2 runs in the cloud as well. Below are 
the input commands to the version of Macaulay2 that runs online. The output is omitted. 

Macaulay2, version 1.6 

11 : loadPackage "PHCpack"; 

12 : help(PHCpack); 

13 : help solveSystem; 

14 : R = CC [x,y,z]; 

15 : S = {x+y+z-1, x"2+y"2, x+y-z-3}; 

16 : solveSystem(S) 
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3 Design and Implementation 


All software in this project is free and open source, as an application of the LAMP stack, where 
LAMP stands for Linux, Apache, MySQL, and Python. Our web server runs Red Hat Linux, 
Apache [6] as the web server, MySQL [9] as the database, and Python as the scripting language. 
Our interest in Python originates in its growing ecosystem for scientific computing [25]. In our 
current implementation we do not take advantage of any web framework. Our web interface 
consists mainly of a collection of Python CGI scripts. 

We distinguish three components in the development of our web interface: the definition of 
the database, the sign up process, and the collection of Python scripts that are invoked as the 
user presses buttons. In the next three paragraphs we briefly describe these three components. 

MySQL is called in Python through the module MySQLdb. The database manages two tables. 
One table stores data about the users, which includes their encrypted passwords and generated 
random names that define the locations of the folders with their data on the server. In the 
other table, each row holds the information about a polynomial system that is solved. In this 
original setup, mathematical data are not stored in the database. Every user has a folder which 
is a generated random 40-character string. With every system there is another generated 40- 
character string. The Python scripts do queries to the database to locate the data that is needed. 

When connecting to the server, the index. html leads directly to the Python script that prints 
the first login screen. The registration script sends an email to the first time user and an activation 
script runs when that first time user clicks on the link received in the registration email. 

The third component of the web interface consists of the Python scripts that interact with 
the main executable program, the phc built with the code in PHCpack. Small systems are solved 
directly. Larger systems are placed in a queue that is served by compute servers. 


4 Solving by Polynomial Homotopy Continuation 

When applying polynomial homotopy continuation methods to solve polynomial systems, we 
distinguish two different approaches. The first is the application of a blackbox solver, and the 
second is a scripting interface. 

4.1 Running a Blackbox Solver 

For the blackbox solver, the polynomials are the only input. The parameters that control the 
execution options are set to work well on a large class of benchmark examples; and/or tuned 
automatically during the solving. While the input is purposely minimal, the output should 
contain various diagnostics and checks. In particular, the user must be warned in case of ill 
conditioning and nearby singularities. The form of the output should enable the user to verify 
(or falsify) the computed results. 

The current blackbox solver phc -b was designed in m for square problems, that is, systems 
with as many equations as unknowns. Polyhedral homotopies |15] m are optimal for sparse 
polynomial systems. This means that the mixed volume is a sharp root count for generic problems 
(phc -b calls MixedVol [12] for a fast mixed volume computation). Every path in a polyhedral 
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homotopy ends at an isolated root, except for systems that have special initial forms [16]. For a 
survey, see e.g. GDI- 

Special cases are polynomials in one variable, linear systems, and binomial systems. A bino¬ 
mial system has exactly two monomials in every polynomial. Its isolated solutions are determined 
by a Hermite normal form of the exponent vectors. Its positive dimensional solution sets are 
monomial maps [1]. 

A more general blackbox solver should operate without any assumptions on the dimension 
of the solution sets. Inspired by tropical algebraic geometry, and in particular its fundamental 
theorem m, we can generalize polyhedral homotopies for positive dimensional solution sets, as 
was done for the cyclic n-roots problems in 0. This general polyhedral method computes 
tropisms based on the initial forms and then develops Puiseux series starting at the solutions of 
the initial forms. 

4.2 The Scripting Interface phcpy 

The Python package phcpy replaces the input and output files by persistent objects. Instead 
of the command line interface of phc with its interactive menus, the user of phcpy runs Python 
scripts, or calls functions from the phcpy modules in an interactive Python shell. Version 0.1.4 
of phcpy is described in [32], 

The current version 0.2.5 exports most of the tools needed to compute a numerical irreducible 
decomposition [29]. With a numerical irreducible decomposition one gets all solutions, the isolated 
solutions as well as the positive dimensional sets. The latter come in the form of as many generic 
points as the degree of each irreducible component, satisfying as many random hyperplanes as 
the dimension of the component. 


5 Pattern Matching with a Database 

Solving a system of polynomial equations for the first time with polynomial homotopy continua¬ 
tion happens in two stages. In the first stage, we construct and solve a simpler system than the 
original problem. This simpler system serves as a start system to solve the original system in the 
second stage. Polynomial homotopy continuation methods deform systems with known solutions 
into systems that need to be solved. Numerical predictor-corrector methods track solution paths 
from one system to another. 

In many applications polynomial systems have natural parameters and often users will present 
systems with the same input patterns. 

5.1 The Classification Problem 

If we solve polynomial systems by homotopy continuation, we first solve a similar system, a 
system with the same monomial structure, but with generic coefficients. If we could recognize the 
structure of a new polynomial system we have to solve, then we could save on the total solving 
time of the new polynomial system, because we skip the first stage of solving the start system. 
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Furthermore, we could give a prediction of the time it will take to solve the new system with the 
specific coefficients. 

Giving a name such as cyclic n-roots as a query to a search engine is likely to come up 
with useful results because this problem has been widely used to benchmark polynomial system 
solvers. But if one has only a particular formulation of a polynomial system, then one would like 
to relate the particular polynomials to the collection of polynomial systems that have already 
been solved. Instead of named polynomial systems, we work with anonymous mathematical data 
where even the naming of the variables is not canonical. 

For example, the systems 

f x 2 + xy 2 — 3 = 0 , f 3 + 2ab 2 = 0 . 

( 2x 2 » + 5 = 0 “ d \6 2 —5 + 2cA = 0 (1) 

must be considered isomorphic to each other. 

The support set of a polynomial is the set of all exponent tuples of the monomials that appear 
in the polynomial with nonzero coefficient. For the systems in the specific example above, the 
sets of support sets are 


{{( 2 , 0 ),( 1 , 2 ),( 0 , 0 )},{( 2 , 1 ),( 0 , 0 )}} ( 2 ) 

and {{(0,0), (1,2)}, {(0,2), (0,0), (2,1)}}. (3) 

Definition 5.1. We say that two sets of support sets are isomorphic if there exists a permutation 
of their equations and variables so that the sets of support sets are identical. 

The problem is then to determine whether the sets of support sets of two polynomial sys¬ 
tems are isomorphic to each other. This problem is related to the isomorphism of polynomials 
problem [23]. Algorithms in multivariate cryptology apply Grobner basis algorithms [10| and 
graph-theoretic algorithms [5]. 

5.2 The Graph Isomorphism Problem 

If we encode the sets of support sets of a polynomial system into a graph, then our problem reduces 
to the graph isomorphism problem, for which practical solutions are available m and accessible 
in software m- The problem of determining whether two sets of support sets are isomorphic 
is surprisingly nontrivial. We begin with some theoretical considerations before moving on to 
implement at ion. 

Definition 5.2. The graph isomorphism problem asks whether for two undirected graphs F, G 
there is a bijection <f> between their vertices that preserves incidence-i.e. if a and b are vertices 
connected by an edge in F, then </>(a) and <f>(b) are connected by an edge in G. 

Proposition 5.3. The problem of determining whether two sets of support sets are isomorphic 
is equivalent to the graph isomorphism problem. 

Proof. We will give a constructive embedding in both directions. 

(D) We start by showing that graph isomorphism can be embedded in checking isomorphism 
of sets of support sets. Recall that the incidence matrix of a graph G = ( V, E) is a matrix with 
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j^V rows and H^E columns where the ( i,j) th entry is 1 if the i th vertex and j th edge are incident 
and 0 otherwise. It is straightforward to show that two graphs are isomorphic if and only if their 
incidence matrices can be made equal by rearranging their rows and columns. 

Now suppose we have a graph with incidence matrix A = (ay). Construct a polynomial by 
considering the rows as variables and the columns as monomials. To be precise, set 

pa = Y,TI x ?‘ 

j i 

For example, if a graph has incidence matrix 

1 M 

10, then pa = x\X 2 + x\Xj, 

0 1 ) 

Switching two columns of the matrix corresponds to reordering the monomials in the sum; switch¬ 
ing two rows corresponds to a permutation of variables. Therefore, to determine if two graphs 
are isomorphic, one may find their incidence matrices, form polynomials from them, and check 
if the polynomials (thought of as single-polynomial systems) are the same up to permutation of 
variables. Finally, it is important to note that this is a polynomial time reduction. To create the 
polynomial, start at the first column and iterate through all of the columns. When done with the 
first column, we have created the first monomial. Repeat this process for every column, touching 
every entry in the matrix exactly once. Given an incidence matrix of a graph G = (V, E ) with 
jj^V = n, #E = m, converting to a polynomial requires 0{nm) operations. 



(5) 


(4) 


equations 


monomials 


exponents 


variables 



Figure 1: Basic graph from system 

(C) We now show the other direction: checking isomorphism of sets of support sets can be 
embedded in graph isomorphism. To do so, we give a way of setting up a graph from a set of 
support sets, which is most easily seen via an example. The graph for the system {x 2 + xy, x 3 y 2 + 
x + 1} is shown in Figure [TJ As can be seen from the diagram, we build a graph beginning with 
a root node to ensure the graph is connected. We attach one node for each equation, then attach 
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nodes to each equation for its monomials. Finally we put a row of variable nodes at the bottom 
and connect each monomial to the variables it contains via as many segments as the degree of 
the variable in the monomial. 

At this point our graph setup is inadequate, as two different systems could have isomorphic 
graphs. Consider the graphs for {x, y} and {a; 4 }, seen in Figure [2j Even though the systems 
are different, the graphs are clearly isomorphic. To remedy this, self-loops are added in order 
to partition the graph into root, equations, monomials, and variables by putting one self-loop at 
the root node, two on each equation node, etc. This graph will uniquely represent our system 
since any automorphism will be a permutation of nodes within their partitions. Partitioning by 
self-loops is possible because there are no self-loops in the initial setup. The graph in Figure Q] is 
drawn without the self-loops for the sake of readability. 


equations 

monomials 

exponents 

variables 



f root 


x 4 ° 


it 

it 

n 
n 
X • 


Figure 2: Isomorphic graphs 

The time to complete this process is easily seen to be polynomial in the number of monomials 
and the sum of their total degrees. 

Because we have shown that two sets of support sets are isomorphic if and only if two uniquely 
dependent graphs are isomorphic, and had previously shown that two graphs are isomorphic if 
and only if two uniquely dependent polynomials are equal up to rearranging variables, we are 
done. Note that this puts our problem in GI, the set of problems with polynomial time reduction 
to the graph isomorphism problem. □ 


5.3 Computing Canonical Graph Labelings With Nauty 

The graph isomorphism problem is one of the few problems which is not known to be P or NP- 
complete m- Although there is no known worst-case polynomial time algorithm, solutions which 
are fast in practice exist. One such solution is nauty [21], a program which is able to compute 
both graph isomorphism and canonical labelings. More information on nauty and other software 
for solving graph isomorphism can be found in [22] , 

Because nauty can incorporate the added information of ordered partitions, we chose to revise 
our setup to take advantage of this and minimize the number of nodes. Instead of using self¬ 
loops to partition the graph, we specify equations, monomials, and variables to be three of the 
partitions. We then check for which exponents occur in the system, and attach nodes representing 
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equations 


monomials 


exponents 

variables 



Figure 3: Nauty graph setup 


these exponents to the variables where appropriate. Instead of using a sequence of nodes and 
edges to record the degree as we did in the proof, we instead attach monomials to these exponent 
nodes. We group these exponent nodes so that all the nodes representing the lowest exponent 
are in the first partition, all representing the second lowest are in the second partition, etc. The 
setup is shown in FigureEl where the different shapes of the exponent nodes represent the distinct 
partitions. 

If we couple this graph with an ordered list of the exponents that occur in the system, we once 
again have an object that uniquely represents the isomorphism class of our system. In addition, 
nauty works by computing the automorphisms of a graph and determining one of them to be 
canonical; the automorphisms must respect the partition, hence partitioning reduces the number 
of possibilities to check. So not only does this use fewer vertices and edges than our previous 
setup, but the use of partitions speeds up nauty’s calculations. 

Another advantage to using nauty is that it computes generators of the automorphism group 
of the graph. Some of these generators permute the equations and monomials. These are unim¬ 
portant and may be discarded, but the generators that permute the variables are quite useful, 
as PHCpack has the ability to take advantage of symmetry in polynomial systems. If a system 
remains the same under certain changes of variables, runtime may be significantly decreased if 
this symmetry is passed to phc. 

It is worth noting that by Frucht’s theorem m and our Proposition 15.31 we immediately 
obtain that for any group G there is a polynomial system with G as its symmetry structure. If 
we want to actually find systems with particular structures this is fairly impractical, however, as 
the proof of Frucht’s theorem uses the Cayley graph which has a node for every group element, 
meaning that a system built with this method would have a variable for every group element. 

We are primarily interested in using nauty in the context of storing tuples of support sets in a 
database. Because of this, the canonical labeling feature is much more useful than the ability to 
check graph isomorphism-looking up a system in a database ought not be done by going through 
a list of systems and querying nauty as to whether they are isomorphic to the one we are looking 
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up, since this would be highly inefficient. Instead we simply parse the system into its graph form 
(including partition data), pass this information to nauty and compute the canonical labeling, 
and attach the exponent data, as a string, to nauty’s output. This process gives us a string 
that uniquely corresponds to the isomorphism class of the system, which we can then store in a 
database. 

5.4 Benchmarking the Canonization 

Timings reported in this section were done on a 3.5 GHz Intel Core i5 processor in a iMac 
Retina 5K running version 10.10.2 of Mac OS X, with 16 GB RAM. Scripts are available at 
https://github.com/sommars/PolyGraph. 

In order to design an intelligent storage system, it is necessary to know an upper bound of the 
length of the string. As polynomial systems can be arbitrarily large, leading to arbitrarily long 
strings from nauty, we chose to analyze a number of well known systems to act as benchmarks, 
giving us an idea of an upper bound. First, consider the cyclic n-root polynomial systems. These 
systems consist of n equations: Y^i= 1 1~I)=i x j mod n = 0 for k from 0 to n — 2, and x\ ... x n — 1 = 0. 
For example, the cyclic-3 system is 


X\ + x 2 + x 3 = 0 

x\X2 + X2X3 + X3X1 = 0 ( 6 ) 

X1X2X3 — 1 = 0 . 


These systems have a lot of symmetry, so they are an interesting case for us to benchmark this 
process. We tested against both small and large values of n to gain a thorough understanding of 
how this system will be canonized. Experiments in small dimensions are summarized in Table [TJ 
The computation time, an average of three trials, is quite fast, though the length of the string 
grows quickly. In comparison, note that the calculation of the root counts in the blackbox solver 
of PHCpack (which includes the mixed volume computation) for the cyclic 10-roots problems 
takes 48.8 seconds. 


n 

time 

#nodes 

^characters 

4 

0.006 

29 

526 

6 

0.006 

53 

1,256 

8 

0.006 

85 

2,545 

10 

0.007 

125 

5,121 

12 

0.007 

173 

8,761 


Table 1: Cyclic n-root benchmark for small n. For each n, its takes only milliseconds. We list 
the number of nodes in the graph and the length of the canonical form. 

For large instances of cyclic n-roots, the exponential growth of the number of solutions in¬ 
creases so much that we may no longer hope to be capable of computing all solutions. Nevertheless, 
with GPU acceleration [?], we can manage to track a limited number of paths. Table [3] illustrates 
the relationship between the dimension n, the time, and the sizes of the data. 
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n 

time 

#nodes 

^characters 

16 

0.010 

293 

20,029 

32 

0.045 

1,093 

168,622 

48 

0.265 

2,405 

601,702 

64 

1.200 

4,229 

1,427,890 

80 

4.316 

6,565 

2,778,546 

96 

15.274 

9,413 

4,784,390 

112 

38.747 

12,773 

8,595,408 

128 

80.700 

16,645 

13,094,752 


Table 2: For larger values of the dimension of the cyclic n-root problem, times and sizes of the 
data start to grow exponentially. 

Another interesting class of benchmark polynomial systems that we can formulate for any 
dimension is the computation of Nash equilibria [8]. We use the formulation of this problem as 
in [23]. The number of isolated solutions also grows exponentially in n. Table [3] summarizes our 
computational experiments. As before, the times represent an average of three trials. 


n 

time 

#nodes 

^characters 

4 

0.006 

47 

977 

5 

0.006 

98 

2,325 

6 

0.007 

213 

7,084 

7 

0.013 

472 

18,398 

8 

0.054 

1,051 

51,180 

9 

0.460 

2,334 

134,568 

10 

4.832 

5,153 

331,456 

11 

73.587 

11,300 

872,893 

12 

740.846 

24,615 

2,150,512 


Table 3: The cost of the canonization for the Nash equilibria polynomial systems for increasing 
dimension n, with the running time expressed in seconds, the number of nodes, and the size of 
the canonical form. 

Compared to the cyclic n-roots problem, the dimensions in Table [3] are just as small as in 
Table [I] but the time and sizes grow much faster for the Nash equilibria than for the cyclic 
n-roots. We suspect two factors. First, while structured, the Nash equilibria systems are not as 
sparse as the cyclic n-roots problems. Second, unlike the cyclic n-roots problem every equation 
in the Nash equilibria system has the same structure, so the full permutation group leaves the 
sets of support sets invariant. 

We end with a system formulated by S. Katsura, see [3J and m- Table [4] shows the cost of 
the canonization of this system. Because there is no symmetry in the support sets, the cost of 
the canonization increases not as fast in the dimension n as with the other two examples. 

The actual cost of the computation of the canonical form may serve as an initial estimate on 
the cost of solving the system. 
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n 

time 

#nodes 

^characters 

25 

0.020 

929 

24,906 

50 

0.090 

3,411 

112,654 

75 

0.546 

7,454 

254,770 

100 

1.806 

13,061 

495,612 

125 

4.641 

20,229 

793,662 

150 

10.860 

28,961 

1,157,498 

175 

21.194 

39,254 

1,587,115 

200 

52.814 

51,111 

2,082,562 

225 

98.118 

64,529 

2,643,891 


Table 4: The cost of the canonization of the Katsura system for increasing dimension n. The 
number of solutions for the n-dimensional version of the system equals 2 n . 

5.5 Storing Labelings in a Database 

There are many different ways to design a database to store a given set of information. We do 
not contend that this is necessarily the best way, but it is certainly an effective way of storing our 
uniquely identifying information that will lead to fast lookup. Consider the explicitly described 
schema in Figure [3J 


Section 

Name 

Datatype 

Description 


njnodejvariable 

INT 

Number of variables nodes 

Graph Nodes 

ri-nodejmonomial 

INT 

Number of monomial nodes 

runode-equation 

INT 

Number of equation nodes 


njnode-degree 

INT 

Number of degree nodes 

Degree Set 

n-degree 

INT 

Sum of all degrees 

degrees 

VARCHAR 

Set of all degrees 


graph-length 

INT 

Length of complete graph 

Graph 

graph-filename 

VARCHAR 

Complete graph in file 

Polynomial Info 

poly-filename 

VARCHAR 

Information and Reference file for poly¬ 
nomial system 


Figure 4: Database structure for polynomial system graph 

Each of the data elements in it are used to partition the database from the previous elements, 
so it can be seen as a B-tree structure. For identifying whether or not a set of support sets is 
already in our database, this would lead to a search time of O(logn) [7|. 


6 Conclusions 

The practical considerations of offering a cloud service to solve polynomial systems with polyno¬ 
mial homotopy continuation led to the classification problem of polynomial systems. To solve this 
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problem, we linked the isomorphism problem for sets of support sets to the graph isomorphism 
problem and applied the software nauty. 

Although no polynomial time algorithm are known to solve the graph isomorphism problem, 
we presented empirical results from benchmark polynomial systems that the computation of a 
canonical form costs much less than solving a polynomial system. 
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